#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>

//#define G 2.959572553206070e-004
#define G 4*3.141592653589793*3.141592653589793
#define diretorioWindows "C:\\dados\\"

void jacobiana3C(double* j, double m1,double m2,double m3,double r1x,double r1y,double r2x,double r2y,double r3x,double r3y,int tamanho);
void jacobianaK(double* j, double m1,double m2,double rx,double ry,int tamanho);
void abrirArquivos(char SO,char* metodo,double tFinal,double n);
double* carregarArquivosLinha(int dimensao,int tamanho);
double* F(double *J,double* delta,int tamanho);
double* produtoEscalarVetor(double* v,double a,int tamanho);
double* somaVetores(double* v,double* v2,int tamanho);
void imprimir(double* v,int tamanho);
double calculaNorma(double *v,int dimensao);

FILE *f[2];

int main(void)
{
	double tFinal,n,*z,*J,*R1,*R2,*R3,*R4,*R5,*R6,dt,*delta,*temp,*norma,*y,*Y;
	int i=0,dimensao,corpos,T,tamanho;
	char aux[100];

	tFinal=1000;
	n=18;
	dimensao=2;
	corpos=1;
	tamanho=2*corpos*dimensao;
	T=pow(2,n)+1;
	//T=2;
	dt=tFinal/pow(2,n);

	abrirArquivos('w',"kepler",tFinal,n);
	J=(double *)calloc(pow(tamanho,2), sizeof(double));
	delta=(double *)calloc(tamanho, sizeof(double));
	delta=(double *)calloc(tamanho, sizeof(double));
	temp=(double *)calloc(tamanho, sizeof(double));
	norma=(double*)calloc(T,sizeof(double));
	y=(double*)calloc(T,sizeof(double));
	Y=(double*)calloc(T,sizeof(double));

	sprintf(aux,"%0.15lf ",Y[0]);
	fputs(aux,f[2]);
	fputs("\n",f[2]);

	//delta inicial aleatorio
	srand( (unsigned)time(NULL) );
	for(i=0;i<tamanho;i++)
	{
		delta[i]=rand()/(double)RAND_MAX;
		delta[i]=delta[i]-0.5;
	}
	//delta[0]=0.1;
	//delta[3]=0.1;

	// norma delta inicial nao normalizado
	norma[0]=calculaNorma(delta,tamanho);
	//normalizacao delta inicial
	for(i=0;i<tamanho;i++)
	{
		delta[i]=delta[i]/norma[0];
	}

	//norma delta inicial normalizado
	norma[0]=calculaNorma(delta,tamanho);

	puts("");

	for(i=1;i<T;i++)
	{
		z=carregarArquivosLinha(dimensao,tamanho);
		jacobianaK(J,1,3.003167261575587e-006,z[0],z[1],tamanho);

		//resolver delta por RK45
		R1=produtoEscalarVetor(F(J,delta,tamanho),dt,tamanho);
		temp=somaVetores(produtoEscalarVetor(R1,0.25,tamanho),delta,tamanho);

		R2=produtoEscalarVetor(F(J,temp,tamanho),dt,tamanho);
		temp=somaVetores(produtoEscalarVetor(R1,3.0/32,tamanho),delta,tamanho);
		temp=somaVetores(produtoEscalarVetor(R2,9.0/32,tamanho),temp,tamanho);

		R3=produtoEscalarVetor(F(J,temp,tamanho),dt,tamanho);
		temp=somaVetores(produtoEscalarVetor(R1,1932.0/2197,tamanho),delta,tamanho);
		temp=somaVetores(produtoEscalarVetor(R2,-7200.0/2197,tamanho),temp,tamanho);
		temp=somaVetores(produtoEscalarVetor(R3,7296.0/2197,tamanho),temp,tamanho);

		R4=produtoEscalarVetor(F(J,temp,tamanho),dt,tamanho);
		temp=somaVetores(produtoEscalarVetor(R1,439.0/216,tamanho),delta,tamanho);
		temp=somaVetores(produtoEscalarVetor(R2,-8.0,tamanho),temp,12);
		temp=somaVetores(produtoEscalarVetor(R3,3680.0/513,tamanho),temp,12);
		temp=somaVetores(produtoEscalarVetor(R4,-845.0/4104,tamanho),temp,12);

		R5=produtoEscalarVetor(F(J,temp,tamanho),dt,tamanho);
		temp=somaVetores(produtoEscalarVetor(R1,-8.0/27,tamanho),delta,12);
		temp=somaVetores(produtoEscalarVetor(R2,2.0,tamanho),temp,12);
		temp=somaVetores(produtoEscalarVetor(R3,-3544.0/2565,tamanho),temp,12);
		temp=somaVetores(produtoEscalarVetor(R4,1859.0/4104,tamanho),temp,12);
		temp=somaVetores(produtoEscalarVetor(R5,-11.0/40,tamanho),temp,12);

		R6=produtoEscalarVetor(F(J,temp,tamanho),dt,tamanho);
		temp=produtoEscalarVetor(R1,16.0/135,tamanho);
		temp=somaVetores(produtoEscalarVetor(R3,6656.0/12825,tamanho),temp,12);
		temp=somaVetores(produtoEscalarVetor(R4,28561.0/56430,tamanho),temp,12);
		temp=somaVetores(produtoEscalarVetor(R5,-9.0/50,tamanho),temp,12);
		temp=somaVetores(produtoEscalarVetor(R6,2.0/55,tamanho),temp,12);

		delta=somaVetores(delta,temp,tamanho);
		norma[i]=calculaNorma(delta,tamanho);

		y[i]=((i-1)/i)*y[i-1]+2*log(norma[i]/norma[i-1]);
		Y[i]=((i-1)*Y[i-1]+y[i])/i;
		//y[i+1]=(i/(i+1))*y[i]+2*log(norma[i+1]/norma[i]);
		//Y[i+1]=((i)*Y[i]+y[i+1])/(i+1);

		sprintf(aux,"%0.15lf ",Y[i]);
		fputs(aux,f[2]);
		fputs("\n",f[2]);
		//imprimir(delta,12);
	}
	imprimir(y,10);

	fclose(f[0]);
	fclose(f[1]);
	fclose(f[2]);

	return 0;
}

void imprimir(double* v,int tamanho)
{
	int i;
	for(i=0;i<tamanho;i++)
	{
		printf("%0.10lf\n",v[i]);
	}
	puts("");
}

double* somaVetores(double* v,double* v2,int tamanho)
{
	int i;
	double *v3=calloc(tamanho,sizeof(double));

	for(i=0;i<tamanho;i++)
	{
		v3[i]=v[i]+v2[i];
	}
	return v3;
}

double* produtoEscalarVetor(double* v,double a,int tamanho)
{
	int i;
	double* v2=0;

	v2=calloc(tamanho,sizeof(double));
	for(i=0;i<tamanho;i++)
	{
		v2[i]=a*v[i];
	}
	return v2;
}

double* F(double *J,double* delta,int tamanho)
{
	int i,j;
	double* F=(double *)calloc(tamanho, sizeof(double)),temp;
	for(i=0;i<tamanho;i++)
	{
		temp=J[i*tamanho]*delta[0];
		for(j=1;j<tamanho;j++)
		{
			temp+=J[i*tamanho+j]*delta[j];
		}
		F[i]=temp;
	}
	return F;
}

void abrirArquivos(char SO,char* metodo,double tFinal,double n)
{
	char *barra,*diretorio,*temp,*raiz;

	raiz=(char*)calloc(1000,sizeof(char));
	diretorio=(char*)calloc(1000,sizeof(char));
	temp=(char*)calloc(1000,sizeof(char));

	if(SO=='w')
	{
		barra="\\";
	}
	else
	{
		barra="/";
	}

	strcpy(raiz,diretorioWindows);
	strcat(raiz,metodo);
	strcat(raiz,barra);
	sprintf(temp, "%0.0lf", tFinal);
	strcat(raiz,temp);
	strcat(raiz,barra);
	sprintf(temp, "%0.0lf", n);
	strcat(raiz,temp);
	strcat(raiz,barra);

	strcpy(diretorio,raiz);
	strcat(diretorio,"r.txt");
	f[0]=fopen(diretorio,"r");

	strcpy(diretorio,raiz);
	strcat(diretorio,"p.txt");
	f[1]=fopen(diretorio,"r");

	strcpy(diretorio,raiz);
	strcat(diretorio,"megno.txt");
	f[2]=fopen(diretorio,"w");

	//strcpy(diretorio,raiz);
	//strcat(diretorio,"r1.txt");
	//f[0]=fopen(diretorio,"r");

	//strcpy(diretorio,raiz);
	//strcat(diretorio,"r2.txt");
	//f[1]=fopen(diretorio,"r");

	//strcpy(diretorio,raiz);
	//strcat(diretorio,"r3.txt");
	//f[2]=fopen(diretorio,"r");

	//strcpy(diretorio,raiz);
	//strcat(diretorio,"p1.txt");
	//f[3]=fopen(diretorio,"r");

	//strcpy(diretorio,raiz);
	//strcat(diretorio,"p2.txt");
	//f[4]=fopen(diretorio,"r");

	//strcpy(diretorio,raiz);
	//strcat(diretorio,"p3.txt");
	//f[5]=fopen(diretorio,"r");
}

double* carregarArquivosLinha(int dimensao,int tamanho)
{
	int k=0;
	char *temp,*temp2;
	double* z;

	z=(double*)calloc(tamanho,sizeof(double));
	temp=(char*)calloc(1000,sizeof(char));
	temp2=(char*)calloc(1000,sizeof(char));

	for(k=0;k<tamanho/2;k++)
	{
		fgets(temp,30,f[k]);
		int i=0,j=0;
		do{
			temp2[i]=temp[i];
			i++;
		}while(temp[i]!=' ');
		z[2*k]=atof(temp2);
		do{
			temp2[j]=temp[i];
			i++;
			j++;
		}while(temp[i]!=' ');
		z[2*k+1]=atof(temp2);
	}
	return z;
}

void jacobianaK(double* j, double m1,double m2,double rx,double ry,int tamanho)
{
	double G2=(double)G;
	double t2 = 1.0/m2;
	double t3 = rx*rx;
	double t4 = ry*ry;
	double t5 = t3+t4;
	double t6 = 1.0/(pow(t5,(5.0/2.0)));
	double t7 = G2*m1*m2*rx*ry*t6*3.0;
	double t8 = 1.0/(pow(t5,(3.0/2.0)));
	//int i,k;

	j[0*tamanho+2]=t2;
	j[1*tamanho+3]=t2;
	j[2*tamanho+0]=-G2*m1*m2*t8+G*m1*m2*t3*t6*3.0;
	j[2*tamanho+1]=t7;
	j[3*tamanho+0]=t7;
	j[3*tamanho+1]=-G2*m1*m2*t8+G*m1*m2*t4*t6*3.0;

	/*for(i=0;i<4;i++)
	{
		for(k=0;k<4;k++)
		{
			printf("%3.6lf ",j[i*4+k]);
		}
		puts("");
	}
	puts("");*/
}

void jacobiana3C(double* j, double m1,double m2,double m3,double r1x,double r1y,double r2x,double r2y,double r3x,double r3y,int tamanho)
{
	double G2=(double)G;
	double t70 = r1x-r2x;
	double t71 = r1y-r2y;
	double t72 = r1x-r3x;
	double t73 = r1y-r3y;
	double t78 = r1x*2;
	double t83 = r2x*2;
	double t74 = t78-t83;
	double t75 = pow(t70,2);
	double t76 = pow(t71,2);
	double t77 = t75+t76;
	double t85 = r3x*2;
	double t79 = -1*t78+t85;
	double t80 = pow(t72,2);
	double t81 = pow(t73,2);
	double t82 = t80+t81;
	double t84 = pow((1/t77),(5.0/2.0));
	double t86 = r1y*2;
	double t87 = pow((1/t82),(5.0/2.0));
	double t88 = pow((1/t77),(3.0/2.0));
	double t89 = pow(t74,2);
	double t90 = G2 * m1 * m2 * t84 * t89* 3.0/4.0;
	double t91 = r2y*2;
	double t92 = t86-t91;
	double t93 = G2*m1*m2*t74*t84*t92*3.0/4.0;
	double t94 = pow((1.0/t82),(3.0/2.0));
	double t95 = t78-t85;
	double t96 = r3y*2.0;
	double t97 = t86-t96;
	double t98 = G2*m1*m3*t87*t95*t97*(3.0/4.0);
	double t99 = t93+t98;
	double t100 = G2*m1*m2*t88;
	double t101 = G2*m1*m3*t94;
	double t102 = pow(t92,2);
	double t103 = G2*m1*m2*t102*t84*(3.0/4.0);
	double t104 = pow(t97,2);
	double t105 = G2*m1*m3*t104*t87*(3.0/4.0);
	double t106 = t100-t90;
	double t107 = r2x-r3x;
	double t108 = r2y-r3y;
	double t109 = t83-t85;
	double t110 = pow(t107,2);
	double t111 = pow(t108,2);
	double t112 = t110+t111;
	double t113 = pow((1.0/t112),(5.0/2.0));
	double t114 = pow((1.0/t112),(3.0/2.0));
	double t115 = pow(t109,2);
	double t116 = G2*m2*m3*t113*t115*(3.0/4.0);
	double t117 = t91-t96;
	double t118 = G2*m2*m3*t109*t113*t117*(3.0/4.0);
	double t119 = t100-t103;
	double t120 = t118+t93;
	double t121 = G2*m2*m3*t114;
	double t122 = pow(t117,2);
	double t123 = G2*m2*m3*t113*t122*(3.0/4.0);
	double t124 = pow(t95,2);
	double t127 = G2*m1*m3*t124*t87*(3.0/4.0);
	double t125 = t101-t127;
	double t126 = -1*t116+t121;
	double t128 = t101-t105;
	double t129 = t121-t123;
	double t130 = t118+t98;

	int colunas=3*2*2;

	j[0*colunas+6]=1/m1;
	j[1*colunas+7]=1/m1;
	j[2*colunas+8]=1/m2;
	j[3*colunas+9]=1/m2;
	j[4*colunas+10]=1/m3;
	j[5*colunas+11]=1/m3;
	j[6*colunas+0]=t90-G2*m1*m2*t88-G2*m1*m3*t94+G2*m1*m3*pow(t79,2)*t87*(3.0/4.0);
	j[6*colunas+1]=t99;
	j[6*colunas+2]=t106;
	j[6*colunas+3]=-t93;
	j[6*colunas+4]=t125;
	j[6*colunas+5]=-t98;
	j[7*colunas+0]=j[0*colunas+1];
	j[7*colunas+1]=-t100-t101+t103+t105;
	j[7*colunas+2]=-t93;
	j[7*colunas+3]=t119;
	j[7*colunas+4]=-t98;
	j[7*colunas+5]=t128;
	j[8*colunas+0]= t106;
	j[8*colunas+1]=-t93;
	j[8*colunas+2]=-t100+t116+t90-G2*m2*m3*t114;
	j[8*colunas+3]=t120;
	j[8*colunas+4]=t126;
	j[8*colunas+5]=-t118;
	j[9*colunas+0]=-t93;
	j[9*colunas+1]=t119;
	j[9*colunas+2]=t120;
	j[9*colunas+3]=-t100+t103-t121+t123;
	j[9*colunas+4]=-t118;
	j[9*colunas+5]=t129;
	j[10*colunas+0]=t125;
	j[10*colunas+1]=-t98;
	j[10*colunas+2]=t126;
	j[10*colunas+3]=-t118;
	j[10*colunas+4]=-t101+t116-t121+t127;
	j[10*colunas+5]=t130;
	j[11*colunas+0]=-t98;
	j[11*colunas+1]=t128;
	j[11*colunas+2]=-t118;
	j[11*colunas+3]=t129;
	j[11*colunas+4]=t130;
	j[11*colunas+5]=-t101+t105-t121+t123;

	//printf("%0.25lf\n",j[3*6+3]);

	/*int i,k;

	for(i=0;i<12;i++)
	{
		for(k=0;k<12;k++)
		{
			printf("%3.6lf ",j[i*colunas+k]);
		}
		puts("");
	}*/
}

double calculaNorma(double *v,int dimensao)
{
	double temp=0;
	int i;

	for(i=0;i<dimensao;i++)
	{
		temp+=v[i]*v[i];
	}
	temp=sqrt(temp);
	return temp;
}



